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Abstract. In these notes we present a pedagogical account of the population dynamics methods 
recently introduced to simulate large deviation functions of dynamical observables in and out of 
equilibrium. After a brief introduction on large deviation functions and their simulations, we review 
the method of Giardina et al. for discrete time processes and that of Lecomte et al. for the continuous 
time counterpart. Last we explain how these methods can be modified to handle static observables 
and extract information about intermediate times. 
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The main achievement of equilibrium statistical mechanics is probably the simplifi- 
cation it offers in the study of static observables in steady state. Indeed, the resolution of 
dynamical equations is replaced by static averages and ensemble approaches. This can 
still be very difficult, but from a conceptual point of view, the problem is much simpler. 
On the other hand, when one is interested in dynamical observables, like currents of par- 
ticles, or in out-of-equilibrium situations, like for glassy or driven systems, such static 
ensemble approaches are not available anymore and there is no general formalism on 
which one can rely. 

For the last ten years, physicists have been interested in large deviation functions 
mainly because they are good candidates to extend the concept of thermodynamic poten- 
tials to out-of-equilibrium situations and to dynamical observables (for a review, see [1]). 
Of course the mere definition of out-of-equilibrium potentials is not useful in itself and 
the challenge is thus to go beyond their construction. To do so, two strategies can be 
followed. First, one can try to derive general properties of large deviation functions. An 
example of success in this direction is provided by the Fluctuation Theorem [2, 3, 4, 5], 
which can be read as a symmetry of large deviation functions and is one of the few 
general results valid out-of-equilibrium. Another strategy is to consider specific exam- 
ples and to compute the large deviation function explicitly. For some simple yet non- 
trivial interacting particle systems, exact computations have been possible (for a review 
see [6]), but one has to rely on numerics for more generic systems. From the algorithmic 
point of view, two paths can be followed. For small system sizes, exact procedures can 
be used (see for instance [8, 9]) but as soon as mesoscopic systems are considered, one 
has to rely on importance sampling approaches. Generalising a procedure followed pre- 
viously to study rare events in chemical reactions [10, 11, 12], Kurchan and co-workers 
developed methods to compute large deviation functions in dynamical systems [13], and 
discrete [14] or continuous [15] time Markov chains. We shall concentrate here on the 
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statistical mechanical aspects and review the methods available for Markov chains. 

Let us consider a system and call {^} its set of configurations. As time goes from to 
a final time t, the system typically jumps into successive configurations "^o, , . . . "^a- at 
distinct times ?i , . . . fj?-. A dynamical observable Q is defined as a sum along the history 
of small contributions <?<;f^_^['4, one for each transition between two successive configu- 
rations. The simplest example of such observable is probably a current of particles Q in 
a one-dimensional lattice gas, which is either incremented of decremented every time a 
particle jumps to the right or to the left, respectively. This contrasts with static observ- 
ables, like the number of particles at a given site, which depend solely on the configura- 
tion of the system at a given time. We will see below that the algorithms used to obtain 
the large deviation functions slightly differ in these two cases. To characterise the fluc- 
tuations of the observable Q, the first thing one can do is to extend the microcanonical 
approach to the space of trajectories and compute the corresponding macrostate entropy 

s{q)^\imhogP[Q{t)^qt]. (1) 
t 

However one knows from usual statistical mechanics that working in the microcanonical 
ensemble is often harder than in the canonical one and we rather introduce a dynamical 
partition function and the corresponding dynamical free energy 

Z(j8,0 = (e-^eW\; ^(^) ^ iimilogZ(/3,0. (2) 

These definitions differ slightly from those used in the mathematical literature, where 
one rather speaks about rate functions —s{q) and cumulants generating functions 
V^(-i8). 

The main purpose of this short review is to explain how one can compute V^(/3) using 
an approach relying on population dynamics. But let us first sketch why direct sampling 
would be inefficient. Consider for simplicity the case where s{q) has a single maximum 
at qo, which satisfies 5(^0) = (for normalisation purpose). Fluctuations around Q — qot 
which occur with probabilities of order one must typically be of order 1 / a/? so that 

P{Q = qt)^ ~ eif(?-9o)'^"(9o) ^ 1 ^ (3) 

whence a probability of larger fluctuations exponentially small in On the other hand, 
the dynamical partition function has a weight e"^^ (with Q ~ qt) exponential in t. As 

a result, there is a competition between the exponentially rare fluctuations and their 
exponential weight such that for j8 of order 1 the trajectories which dominate the average 
in (2) are exponentially rare. A more quantitative way to rephrase this can be read in the 
Legendre relation between s{q) and V^(j8): V^(j8) = max^[s(^) — Pq]. The maximum 
is realised for a value q* which dominates the average in (2). It thus differs from qo 
- which maximises only s{q) - by a factor independent of t. From (3) we see that 
the corresponding trajectories indeed have a probability exponentially small in t. To 
have a good sampling over unbiased simulations, should thus be of order e^, an 
impracticable requirement. Direct sampling is thus a hopeless strategy to observe large 
deviations. 
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We present in this section the method first introduced by Giardina, Kurchan and 
Peliti [14] to simulate cumulant generating functions in discrete time Markov chains. 
In this case, the dynamics is defined by the transition probabilities U{'t^ —>■ = [/■^'^ 
between configurations and the corresponding the master equation reads 

P{^,t)^Y^U.^^>P{^',t-l). (4) 

Conservation of probability enforces the matrix U to be stochastic, i.e. for all 
L"^^"^^' — 1- Starting from a fixed configuration the explicit solution of (4) is 

P{^,t)= £ U^^,_,U^,_,^,_,...U^,^,= [U%^^. (5) 

The dynamical observable Q can be written as a sum over configuration changes Q{t) = 
+ • • • + ^^i-^o- compute the dynamical free energy we first rewrite the 

dynamical partition function as 



(6) 



where we have introduced the matrix [f/jg]-^-^' = U%'%^ie^^^'^"^"' . Let us note that i/a(/3) is 
given by the log of the largest eigenvalue oiU^. A possible strategy, used for instance 
in [7, 8, 9], is thus to compute numerically this eigenvalue. The matrix h however 
exponentially large in the system size, which limits this strategy to small systems. 
The main advantage of this method is to yield a numerical approximation of an exact 
expression - as pointed out in [9] - as opposed to our approach, efficient for large 
systems, but relying on importance sampling. Note that the 'exact' approach can be 
used to check the validity of the importance sampling approach for small system sizes, 
before going to larger ones, as was actually done for the simulations presented in [15]. 

Comparison of expressions (5) and (6) leads one to think that ^f{fi) could be obtained 
from a new dynamics, induced by U^. However, the matrix is not stochastic, as in 
general 

% = L[t^i8]^^' = £C/^^'e-^«^^' ^ 1 (7) 

and we should not understand (6) as a stochastic process with conserved probability but 
as a population dynamics with branching and death, where the population size is not 
constant. To do so, let us define 

UL^,J^^ = ^e-^^-'^'. (8) 

The matrix U' is stochastic and (6) now writes 
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This expression is now closer to (5) as U' is stochastic, and can be interpreted as follows: 
N agents evolve with the stochastic dynamics defined by U' and are replicated with a rate 
when they are in configuration ^. This interpretation is possible as Y<^ depends solely 
on the initial configuration and can thus be interpreted as a configuration-dependent 
reproduction rate. This was less apparent in (6), where factors depend on both 

initial and final configurations. 

This interpretation as a population dynamics can be implemented using a diffusion 
Monte Carlo algorithm. Let us consider an ensemble of A^o agents (A/q ^ 1) evolving in 
the configuration space {'^}. At each time step T ^ T + 1, 

(1) Each agent evolves according to the /3 -modified dynamics Vl^r^i, 

(2) Each agent in configuration ^ is replicated/killed with probability IV, i.e. is re- 
placed by y copies, where 

„_/ L%J + 1 with probability - [y<^ J 

\Ycg\ with probability l-(y^-[y^J) ^'"^ 

Concretely, the agent is replaced by y copies of itself, so that the population size is 
increased by y — 1 (decreased by 1 if j = 0). 

Let us show that the size of the population at time t yields the large deviation function. 
For a given history, the number A^(^, t) of copies in configuration ^ at intermediate 
time T satisfies N{'i^T, t) = U^^c^^_Yc^^_^N{'^-^-i,T — 1), so that for the whole history 

N{%,t) =U!^^^^_Y<^,_,...U!^^^Y^,N{^o,0). (11) 

Consequently, we see from (9) that the population size N{t) = Y^'^^i^it) behaves as 
N{t) /No = ~ e^^^^\ As usual in importance sampling approaches, the ensem- 

ble average (.) has been replaced by an average over a finite number of simulations. 
Whereas in principle correct, this approach is however impracticable since the popu- 
lation size varies exponentially in time and we thus add a third step to the previous 
algorithm: 

(3) After the cloning step, the population is rescaled by a factor to its initial size Nq, 
by uniformly pruning/replicating the agents. 

At each time step T, the rescaling factor is given by = so that Xt...XQ = 

and finally 

V^(j8) = -limilog(Zf...Xo). (12) 

CONTINUOUS TIME DYNAMICS 

For a continuous time dynamics defined by rates ^ 'to") the master equation reads 
a,P(^,?)= £ W(^'^^)P(^',?) -r(^)P(^,?), (13) 
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where r{^) = Y^(^iW{^ is the escape rate from configuration 

There are many different ways of deriving the algorithm presented in the previous 
section and we shall follow here a derivation of the continuous time algorithm slightly 
different from the one we introduced in [15]. The formal solution of (13) reads 

n'^,?) =I^>oL^i...^^_Ji^^if/^ ^^/^-i---/^ ^^1 (14) 

nit \^ t \ n(t \^ f \^-{t-tK)rm W{%^Vi) Wj^K-x^"^) 
p[tK[eK-\,tK-\)---p[ti[^Q,tQ)e ^ r(^o) r[^K-l) ' 

where p{tk\'^k-\ih-\) = r{%_i)exp[-{tk - tk-i)r{'^k-l)] represents the probability 
distribution of the time intervals between jumps. The sum over K corresponds to all 
the possible numbers of jumps between and t, the sum over the ^^'s to the different 
configurations which can be visited. The integrals over account for all the possible 
times at which jumps occur. Last, the ratio ^^'^^'^n'^^-' gives the probability that the 
system goes to configuration when it quits configuration ^k-l- Multiplying the 
second line of (14) by e~^^ yields an explicit formula for Z(P,t): 

Z{P,t)^ LK>oL'^,...'^,_,,'^gdtKjl^dtK-i...j;^dh (15) 
p{tK\'^K-utK-i)---piti\%,to)e-^^-^^y(^^ 

Further introducing the biased rates Wpi"^ ^') = W^tf ^')e~^««"^, the corre- 
sponding escape rates and time distributions between two jumps pp, (15) can be 
rewritten (after some algebra) as 

Z(j8,0 = LK>oL%...'^^_,,<^l!,dtKi;,^dtK-i...lSdti (16) 

Pp{tK\'^K-i,tK-i)---pp{ti\%,to)e-^'~'''^''P^'^^ 



where F(^fc) = e''/3('^*)~''('^*). Z(j8 , t) is thus a weighted sum over all possible trajectories 
generated by the biased rates W^g , where the weights are given by the factors Y 
A first idea which can come to mind is to simply evolve the population with the rates 
Wp, without cloning and to simply average e./o'^'^(''/5('^(^))^''('^(^))) over these trajectories, 
as was proposed in [16]. This however fails as soon as t is large, for the same reason 
as the one described in the introduction: the weight is exponentially large in t and large 
fluctuations of the exponent are exponentially rare. One thus has to use a biased sampling 
to compute the average (16). Following the philosophy of the "Go with the Winner 
methods" [19], the general idea is to stochastically replace a trajectory with weight W 
hy trajectories with weight 1, so that trajectories which high rates are favoured 
whereas those with small weights are not investigated. 

If the re-weighting procedure is made systematic, every time an agent Ca changes of 
configuration at time t", one gets the following algorithm: 

(0) The time is set to f". 
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(1) Ca jumps from its configuration ^ to another configuration with probability 

(2) The time interval A? until the next jump of Ca is chosen from the Poisson law of 
parameter r^(^'). 

(3) The agent Ca is either cloned or pruned with a rate ^C^') = 

a) One computes y = {^') + ej where e is uniformly distributed on [0, 1]. 

b) If y = 0, the copy Ca is erased. 

c) If J > 1, we make y — \ new copies of Cq. 

(4) If J = 0, one agent ^ Ca is chosen at random and copied, while iiy > l,y —I 
agents are chosen uniformly among the + y — 1 agents and erased. We store the 
rescaling factor X = jj^zri- 

To reconstruct the dynamical free energy, we keep track of all X factors 

ylog(Xi...X^) = |log^e-/^2(f)^ ^ _y^(j8) as t^oo (17) 
Once again the step (4) ensures constant population. 

THE CASE OF STATIC OBSERVABLES 

The methods presented above only apply for dynamical observables, which can be de- 
composed as sums of individual contributions over each configuration change. One 
could also be interested in averages of static observables along the histories O = 
/oJto(t) = ~h)o{^k)- In this case, the above procedure simplifies and the 

algorithm is identical apart from two points. 

• First, there is no bias in the rate. The agents are evolved with the unmodified 
Markov rates W(^ ^ <^')- 

• Second, the cloning rate is simply given by 

This can best be seen in formula (15) by replacing the weight e'^'^'^^+i'^k which depends 
on the configurations before and after the jump by which depends solely 

on the configuration '^^^ and can thus be seen as a constant cloning rate for the whole time 
the system spends in the configuration 

INTERMEDIATE TIMES 

As pointed out in [14], the configurations probed along the simulation are representative 
of the typical ones at final time t in the evolution, rather than at intermediate times 
(0 <^ T ^ t). In particular, the weighted average value {0{t))p of a static observable 
0{^{t)) at final time t is obtained in the algorithm by computing the average of O 
among the agents at the end of the simulation. 

In general, the value of {0{T))p at intermediate times <^ t t differs from the 
one at final time, yet {0{z))p is of particular interest since it is representative of 
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configurations visited during most of the weighted evolution (see [20] for examples). 
A way to compute (0(t))^ numerically can be read in its formal expression: 



{0{T))p = l^K>oL^,...'^,_,,<^i:,dtKg'dtK-i...ll'dHO{z) (18) 

One simply runs the same algorithm as before, which generates the bias on trajectories, 

except that whenever an agent arrives at a time t^^ such that t^-i < T < t^, the corre- 
sponding value 0{^k_i) is attached to the agent. Then, each time an agent is cloned, the 
corresponding value of O is copied accordingly. At the end of the simulation, {0{z))p 
is obtained from the average of the values of 0(t) attached to the surviving clones. Of 
course, thanks to the cloning process between T and t, this average differs from the one 
we could have done at the intermediate time r in the simulation. 

The same kind of scheme also applies to compute the weighted average of any 
observable O depending on the whole history of the system and the crucial step is to 
copy the value of the observable when cloning events occur. Note in particular that the 
determination of {0{T))p can be quite noisy since only a few instances of 0{t) have 
survived at time t. In the long time limit one may gain similar information by studying 
i(/o(iTO('^(T)))o [20] which is a less noisy dynamical observable. 



DISCUSSION 



These numerical methods have been applied successfully in many different situations but 
it is important to keep in mind their limitations. First, as pointed out in [17] convergence 
problems are met when the evolution operator is gapless. This can for instance happen 
in systems where the configuration space is unbounded, as in the Zero Range Process 
but will however not be a problem as long as the configuration space remains finite. 




FIGURE 1. Large deviation function i/a(j3)/L (in log scale) for the total current of particles in the 
SSEP at density p = 1/2 (L = 400 sites). The points (+) are the result of the continuous-time numerical 
algorithm. The line is the analytic result (19) vaUd for very large deviations |j3 1 » 1. 
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The other important limitation is the finiteness of the population of agents, which can 
be problematic for very large deviations yielding large cloning rates. For all applications 
considered in [15], we thus ensured that the cloning factor would never grow larger than 
few percent of the total population size and agreement with theory - when at hand - 
was very good. As an example, one can compare (figure 1) for the simple symmetric 
exclusion process (SSEP) the numerical result in the regime of very large deviations 
with the analytical result [18] 

iv^(^) =2cosh^^-2p(l-p)-2^^^^^^ + ^(.-|^l) (19) 

valid for |/3 1 ^ 1, for the total current of particles in the system. Agreement is very good 
although the values of j8 correspond to very large deviations. For large cloning rates, it 
may also be necessary to modify step (4) of the continuous time algorithm so that agents 
are not pruned uniformly but according to the weight they carry since their last change of 
configuration £i'^^k)l''i}('^k)-r{%)] j^yj ^^j. simulations we never ran into this problem. 
In this worst case scenario, the efficiency of the continuous time implementation would 
fall back to that of the discrete time where at every step the whole population is re- 
sampled. 
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